pop_size
## [1] 500
number_of_generations
## [1] 800
gs_vec
## [1]  3  7 13 18 41
fnc_vec
## [1] 1
initial_variance_vec
## [1] 0
initial_knowledge_vec
## [1] 0
average_pdrift_vec
## [1] 0.1 0.2 0.3
variance_pdrift
## [1] 0
time_cost_vec 
## [1] 0.01 0.02 0.03
collective_vec
## [1] 0 1
result_grid

#Static environment with evolving sharing and innovation
load(paste0(name,"/",1,name,".RData"))
final_data_all<-expand.grid(Group_Size=gs_vec,FNC=fnc_vec,average_pdrift=average_pdrift_vec,initial_knowledge=initial_knowledge_vec,initial_variance=initial_variance_vec,time_cost=time_cost_vec,collective=collective_vec,m1=NA,sd1=NA,m2=NA,sd2=NA,m3=NA,sd3=NA,m4=NA,sd4=NA,m5=NA,sd5=NA,m6=NA,sd6=NA,m7=NA,sd7=NA)  
final_data_all[1,]<-final_data[1,]


for (loop_parameters in 1:nrow(final_data)){
 if( file.exists(paste0(name,"/",loop_parameters,name,".RData")) ){
   load(paste0(name,"/",loop_parameters,name,".RData"))
   final_data_all[loop_parameters,]<-final_data[loop_parameters,]
      } else {
     print("File missing")
     }
  
  if (runif(1)<1.2){
  par(mfrow=c(2,4))
  matplot((result[[1]]),type="l",xlab="Generation",ylab="bound_initial_bias",col="gray",ylim=c(lb[1],ub[1]));lines(1:number_of_generations,apply(result[[1]],1,mean),lwd=2)
  matplot((result[[2]]),type="l",xlab="Generation",ylab="bound_theta",col="gray",ylim=c(lb[2],ub[2]));lines(1:number_of_generations,apply(result[[2]],1,mean),lwd=2)
  matplot((result[[3]]),type="l",xlab="Generation",ylab="delta",col="gray",ylim=c(lb[3],ub[3]));lines(1:number_of_generations,apply(result[[3]],1,mean),lwd=2)
  matplot((result[[4]]),type="l",xlab="Generation",ylab="SI use",col="gray",ylim=c(lb[4],ub[4]));lines(1:number_of_generations,apply(result[[4]],1,mean),lwd=2)
  matplot((result[[5]]),type="l",xlab="Generation",ylab="First",col="gray",ylim=c(lb[5],ub[5]));lines(1:number_of_generations,apply(result[[5]],1,mean),lwd=2)
  matplot((result[[6]]),type="l",xlab="Generation",ylab="Perfomance",col="gray",ylim=c(-0.1,1));lines(1:number_of_generations,apply(result[[6]],1,mean),lwd=2)
  matplot((result[[7]]),type="l",xlab="Generation",ylab="Perfomance of first",col="gray",ylim=c(-0.1,1));lines(1:number_of_generations,apply(result[[7]],1,mean),lwd=2)
  Sys.sleep (.2)
  }
}

#Static environment with evolving sharing and innovation
load(paste0(name,"/",1,name,".RData"))
final_data_all<-expand.grid(Group_Size=gs_vec,FNC=fnc_vec,average_pdrift=average_pdrift_vec,initial_knowledge=initial_knowledge_vec,initial_variance=initial_variance_vec,time_cost=time_cost_vec,collective=collective_vec,m1=NA,sd1=NA,m2=NA,sd2=NA,m3=NA,sd3=NA,m4=NA,sd4=NA,m5=NA,sd5=NA,m6=NA,sd6=NA,m7=NA,sd7=NA)  
final_data_all[1,]<-final_data[1,]


for (loop_parameters in 1:nrow(final_data)){
 if( file.exists(paste0(name,"/",loop_parameters,name,".RData")) ){
   load(paste0(name,"/",loop_parameters,name,".RData"))
   final_data_all[loop_parameters,]<-final_data[loop_parameters,]
      } else {
     print("File missing")
     }
  
  if (runif(1)<1.2){
  par(mfrow=c(2,4))
  matplot((result[[1]]),type="l",xlab="Generation",ylab="bound_initial_bias",col="gray",ylim=c(lb[1],ub[1]));lines(1:number_of_generations,apply(result[[1]],1,mean),lwd=2)
  matplot((result[[2]]),type="l",xlab="Generation",ylab="bound_theta",col="gray",ylim=c(lb[2],ub[2]));lines(1:number_of_generations,apply(result[[2]],1,mean),lwd=2)
  matplot((result[[3]]),type="l",xlab="Generation",ylab="delta",col="gray",ylim=c(lb[3],ub[3]));lines(1:number_of_generations,apply(result[[3]],1,mean),lwd=2)
  matplot((result[[4]]),type="l",xlab="Generation",ylab="SI use",col="gray",ylim=c(lb[4],ub[4]));lines(1:number_of_generations,apply(result[[4]],1,mean),lwd=2)
  matplot((result[[5]]),type="l",xlab="Generation",ylab="First",col="gray",ylim=c(lb[5],ub[5]));lines(1:number_of_generations,apply(result[[5]],1,mean),lwd=2)
  matplot((result[[6]]),type="l",xlab="Generation",ylab="Perfomance",col="gray",ylim=c(-0.1,1));lines(1:number_of_generations,apply(result[[6]],1,mean),lwd=2)
  matplot((result[[7]]),type="l",xlab="Generation",ylab="Perfomance of first",col="gray",ylim=c(-0.1,1));lines(1:number_of_generations,apply(result[[7]],1,mean),lwd=2)
  Sys.sleep (.2)
  }
}

print(simulation)
## function (plot_yes, n, initial_variance, initial_knowledge, average_pdrift, 
##     variance_pdrift, initial_bias, theta, delta, false_costs, 
##     time_cost) 
## {
##     p_drift = rnorm(n, average_pdrift, variance_pdrift)
##     base_line = 0.5
##     d_time = rep(NA, n)
##     h = 0.05
##     sigma = 1
##     t = 1
##     predator = runif(1) < base_line
##     if (predator) 
##         intercept = rnorm(n, 0.5 * theta + initial_knowledge + 
##             initial_bias, initial_variance)
##     if (!predator) 
##         intercept = rnorm(n, 0.5 * theta + initial_knowledge - 
##             initial_bias, initial_variance)
##     L = intercept
##     attribute = 1
##     dec = rep(FALSE, n)
##     social_strength = 0
##     new_dec = TRUE
##     if (plot_yes == 1) {
##         L_large = L
##     }
##     while (any(dec == FALSE)) {
##         dec_before <- dec
##         d_time[!(L <= 0 | L >= theta) == dec] = t
##         dec = L <= 0 | L >= theta
##         if (new_dec == TRUE) {
##             if (sum(dec) > 0) {
##                 social_strength = (sum(L[dec] >= theta[dec] * 
##                   0.5) - sum(L[dec] <= theta[dec] * 0.5))
##             }
##         }
##         randomNorm = rnorm(n, 0, 1)
##         evidence = (p_drift + (social_strength * delta)) * h + 
##             (sigma * h^0.5) * randomNorm
##         evidence[dec] = 0
##         L = L + evidence
##         t = t + 1
##         new_dec <- any((dec) != (dec_before))
##         if (plot_yes == 1) {
##             L_large <- cbind(L_large, L)
##         }
##     }
##     d_time[is.na(d_time)] <- t
##     d_time = (d_time - 1) * h
##     final_decision <- L >= theta * 0.5
##     if (predator) 
##         fitness = ifelse(final_decision, 1, -false_costs[2]) - 
##             (d_time * time_cost)
##     if (!predator) 
##         fitness = ifelse(final_decision, 1, -false_costs[1]) - 
##             (d_time * time_cost)
##     sim_data <- data.frame(final_decision, d_time, fitness, predator)
##     if (plot_yes == 1) {
##         matplot(seq(h, t * h, h), t(L_large), type = "l", ylim = c(0, 
##             max(theta)))
##     }
##     return(sim_data)
## }
print(runEA)
## function (numAgents, pop_size, number_of_generations, gif, false_costs, 
##     time_cost, initial_knowledge, initial_variance, average_pdrift, 
##     variance_pdrift, collective) 
## {
##     if (collective == 0) 
##         n_eval = ceiling(pop_size/as.numeric(numAgents) * 25)
##     if (collective == 1) 
##         n_eval = ceiling(pop_size * 10)
##     library(scales)
##     result = NULL
##     init = "dev"
##     genpool = matrix(NA, pop_size, length(ub))
##     if (init == "rand") {
##         for (gene_i in 1:length(ub)) {
##             start_mean_gen <- runif(1, 0.4, 0.6)
##             a = 1
##             b = (a/start_mean_gen) - a
##             chromosom = rbeta(pop_size, a, b)
##             chromosom = rescale(chromosom, c(lb[gene_i], ub[gene_i]), 
##                 c(0, 1))
##             genpool[, gene_i] <- chromosom
##         }
##     }
##     else {
##         inits <- c(0, 1, 0.1, 0, 0)
##         for (gene_i in 1:length(ub)) {
##             genpool[, gene_i] <- inits[gene_i]
##         }
##     }
##     result = rbind(result, c(apply(genpool, 2, mean), 0, 0))
##     fitness_long = NULL
##     results_long = NULL
##     for (gen in 1:(number_of_generations - 1)) {
##         probs = rep(1, pop_size)
##         fitness = rep(0, pop_size)
##         count = rep(0, pop_size)
##         first_perf = rep(NA, n_eval)
##         for (evals in 1:n_eval) {
##             if (collective == 0) {
##                 index = sample(1:pop_size, numAgents, prob = probs)
##             }
##             else {
##                 index = sample(1:pop_size, 1, prob = probs)
##             }
##             count[index] = count[index] + 1
##             probs[index] = probs[index]/5
##             probs <- probs/max(probs)
##             plot_yes = 0
##             deltas = genpool[index, 3]
##             deltas[genpool[index, 4] == 0] = 0
##             deltas[genpool[index, 5] == 1] = 10
##             if (collective == 0) {
##                 sim_result <- simulation(plot_yes, numAgents, 
##                   initial_variance, initial_knowledge, average_pdrift, 
##                   variance_pdrift, initial_bias = genpool[index, 
##                     1], theta = genpool[index, 2], delta = deltas, 
##                   false_costs, time_cost)
##                 fitness[index] <- fitness[index] + sim_result$fitness
##             }
##             else {
##                 sim_result <- simulation(plot_yes, numAgents, 
##                   initial_variance, initial_knowledge, average_pdrift, 
##                   variance_pdrift, initial_bias = rep(genpool[index, 
##                     1], numAgents), theta = rep(genpool[index, 
##                     2], numAgents), delta = rep(deltas, numAgents), 
##                   false_costs, time_cost)
##                 fitness[index] <- fitness[index] + mean(sim_result$fitness)
##             }
##             firsts = sim_result$fitness[sim_result$d_time == 
##                 min(sim_result$d_time)]
##             means = mean(firsts)
##             first_perf[evals] <- firsts[nnet::which.is.max(-abs(firsts - 
##                 means))]
##         }
##         if (gif == 1) {
##             fitness_long = c(fitness_long, (fitness))
##             results_long = rbind(results_long, (genpool))
##         }
##         k = 5
##         winners1 = winners2 = rep(NA, pop_size)
##         for (turns in 1:pop_size) {
##             indexes = sample(1:pop_size, k)
##             winners1[turns] = indexes[which.max(fitness[indexes])]
##         }
##         genpool <- genpool[winners1, ]
##         for (cross in 1:ncol(genpool)) {
##             corssover_rate = 0.01
##             dummy_index <- which(runif(pop_size) < corssover_rate)
##             dummy_index_partner <- sample(pop_size, length(dummy_index))
##             if (length(dummy_index) > 0) {
##                 for (corss_process in 1:length(dummy_index)) {
##                   which_gen <- sample(length(ub), 1)
##                   if (which_gen %in% c(4, 5)) 
##                     which_gen <- c(4, 5)
##                   genpool[c(dummy_index[corss_process], dummy_index_partner[corss_process]), 
##                     which_gen] = genpool[c(dummy_index_partner[corss_process], 
##                     dummy_index[corss_process]), which_gen]
##                 }
##             }
##         }
##         for (mut in 1:ncol(genpool)) {
##             mutation_rate = 0.01
##             dummy_index <- runif(pop_size) < mutation_rate
##             if (mut <= 3) {
##                 genpool[dummy_index, mut] <- genpool[dummy_index, 
##                   mut] + rnorm(sum(dummy_index), 0, (ub - lb)[mut] * 
##                   0.05)
##                 genpool[genpool[, mut] < lb[mut], mut] <- lb[mut] - 
##                   (genpool[genpool[, mut] < lb[mut], mut] - lb[mut])
##                 genpool[genpool[, mut] > ub[mut], mut] <- ub[mut] - 
##                   (genpool[genpool[, mut] > ub[mut], mut] - ub[mut])
##             }
##             else {
##                 genpool[dummy_index, mut] = rbinom(sum(dummy_index), 
##                   1, 0.5)
##                 if (mut == 4) 
##                   genpool[dummy_index, mut + 1][genpool[dummy_index, 
##                     mut] == 1] = 0
##                 if (mut == 5) 
##                   genpool[dummy_index, mut - 1][genpool[dummy_index, 
##                     mut] == 1] = 0
##             }
##         }
##         result = rbind(result, c(apply(genpool, 2, mean), ifelse(collective == 
##             0, mean(fitness)/(n_eval/(pop_size/numAgents)), mean(fitness)/(n_eval/(pop_size))), 
##             mean(first_perf)))
##     }
##     if (gif == 1) {
##         library(gganimate)
##         library(tweenr)
##         library(magick)
##         library(purrr)
##         library(tidyr)
##         library(dplyr)
##         library(cowplot)
##         plot_data = data.frame(ease = rep("cubic-in-out", length(fitness_long)), 
##             fitness_long, x1 = results_long[, 1], x2 = results_long[, 
##                 2], x3 = results_long[, 3], x4 = results_long[, 
##                 4], x5 = results_long[, 5], id = rep(1:pop_size, 
##                 (number_of_generations - 1)), gen = sort(rep(1:(number_of_generations - 
##                 1), pop_size)))
##         number_of_frames = number_of_generations * 2
##         if (number_of_generations > 99) {
##             number_of_frames = round(number_of_generations/2)
##             plot_data = plot_data %>% filter(gen %in% seq(1, 
##                 number_of_generations, 5))
##         }
##         plot_data2 = plot_data %>% gather(type, phen, c(x1, x2, 
##             x3, x4, x5))
##         plot_data2$max = ifelse(plot_data2$type == "x1", ub[1], 
##             ifelse(plot_data2$type == "x2", ub[2], ub[3]))
##         plot_data2$min = ifelse(plot_data2$type == "x1", lb[1], 
##             ifelse(plot_data2$type == "x2", lb[2], lb[3]))
##         plot_data2 = plot_data2 %>% mutate(type = ifelse(type == 
##             "x1", "Bias", ifelse(type == "x2", "Boundary Seperation", 
##             ifelse(type == "x3", "Social Drift", ifelse(type == 
##                 "x4", "SI use", "Follow")))))
##         pnew = ggplot(plot_data2, aes(x = phen, y = fitness_long, 
##             color = as.factor(id))) + geom_point() + transition_states(gen, 
##             1, 1) + facet_wrap(~type, scales = "free") + theme(legend.position = "none") + 
##             labs(title = "Generation: {closest_state}", x = "Phenotype", 
##                 y = "Fitness") + ylim(0, max(plot_data2$fitness_long))
##         anim_save(paste("Output/tournement", "test", ".gif"), 
##             animation = animate(pnew, nframes = number_of_frames, 
##                 width = 1100, height = 600))
##     }
##     cbind(result)
## }